######################################################
# Code to create figures/analyses in Online Appendix D
######################################################

############################
#Packages
############################
library(survey)
library(srvyr)
library(directlabels)
library(rio)
library(tidyverse)
library(ggpubr)

############################
#Import Data and Clean
############################

#######Data
combdata <- import("./Data/pew_cnn_comb.Rda")

# combdata$sdate <- as.Date(combdata$sdate, 
#                           "%Y-%m-%d")       


                      ########################
                      # Figure OE1: % DK
                      ########################

combdata <- combdata %>%
  mutate(dks = case_when(
    rep_dk == 0 & dem_dk == 0 ~ 1, 
    rep_dk == 1 & dem_dk == 0 ~ 2, 
    rep_dk == 0 & dem_dk == 1 ~ 2, 
    rep_dk == 1 & dem_dk == 1 ~ 3), 
    textreme_nodk = case_when(
      too_extreme1 == 1 & dks == 1 ~ 1, 
      too_extreme1 == 2 & dks == 1 ~ 2, 
      too_extreme1 == 3 & dks == 1 ~ 3))
      

#######Declare Survey Design
comb_design <- combdata %>%
  filter(!is.na(weight)) %>%  #This removes some observations with missing weight data
  srvyr::as_survey_design(ids = 1, weights=weight)

########Get weighted means
dk_means <- comb_design %>%
  mutate(dks = factor(dks, 
                      levels=c(1,2,3), 
                      labels=c("No DK", 
                               "1 DK", 
                               "Both DK"))) %>%
  filter(!is.na(dks)) %>%
  group_by(sdate, dks) %>%
  summarize(proportion = srvyr::survey_mean(na.rm=T)) %>%
  mutate(lower = proportion - (1.96*proportion_se), 
         upper = proportion + (1.96*proportion_se), 
         sdate = as.Date(sdate))

comb_design %>%
  mutate(dks = factor(dks, 
                      levels=c(1,2,3), 
                      labels=c("No DK", 
                               "1 DK", 
                               "Both DK"))) %>%
  filter(!is.na(dks)) %>%
  group_by(sponsor, dks) %>%
  summarize(proportion = srvyr::survey_mean(na.rm=T)) 

#######Figure
f1 <- dk_means %>%
  mutate(label =if_else(sdate == max(sdate), as.character(dks), NA_character_)) %>%
  ggplot(aes(x=sdate, y=proportion, shape=dks, linetype=dks)) + 
  geom_point() + 
  geom_smooth(method=lm, se=F, color="black") + 
  scale_y_continuous(limits=c(0,1)) + 
  theme_light() + 
  labs(y = "Proportion Giving Response", 
       x = "Time",
       shape = "Response", 
       linetype = "Response") + 
  scale_x_date(date_breaks= "2 year", date_labels="%Y") + 
  theme(legend.position="top",
        legend.direction="horizontal", 
        legend.background = element_rect()) + 
  scale_shape_manual(values=c(17, 1, 0)) + 
  scale_linetype_manual(values=c(1, 3, 2)) + 
  theme(legend.position=c(0.30,0.95),
        legend.direction="horizontal", 
        legend.background = element_rect(), 
        legend.box.background = element_rect(colour = "black")) 

f1a <- lemon::reposition_legend(f1, 'top left', offset= 0.002 )  

ggsave(plot = f1a, "figure_oe1.png", 
       height=6, width=10, dpi=1200)

                      
################################################################################################
                # Figure OE2: Replication of Figure 1
###############################################################################################

overall_means <- comb_design %>%
  mutate(too_extreme1 = factor(too_extreme1, 
                               levels=c(1,2,3), 
                               labels=c("Neither", "One Party", "Both Parties"))) %>%
  filter(!is.na(too_extreme1)) %>%
  group_by(sdate, too_extreme1) %>%
  summarize(proportion = srvyr::survey_mean(na.rm=T)) %>%
  mutate(lower = proportion - (1.96*proportion_se), 
         upper = proportion + (1.96*proportion_se), 
         sdate = as.Date(sdate), 
         analysis = "With DKs") %>%
  rename(dv = too_extreme1)

overall_nodk <- comb_design %>%
  mutate(textreme_nodk = factor(textreme_nodk, 
                               levels=c(1,2,3), 
                               labels=c("Neither", "One Party", "Both Parties"))) %>%
  filter(!is.na(textreme_nodk)) %>%
  group_by(sdate, textreme_nodk) %>%
  summarize(proportion = srvyr::survey_mean(na.rm=T)) %>%
  mutate(lower = proportion - (1.96*proportion_se), 
         upper = proportion + (1.96*proportion_se), 
         sdate = as.Date(sdate), 
         analysis = "No DKs") %>%
  rename(dv = textreme_nodk)

overall_comb <- bind_rows(overall_means, overall_nodk)

#######Figure

f2 <- overall_comb %>%
  mutate(label =if_else(sdate == max(sdate), as.character(dv), NA_character_)) %>%
  ggplot(aes(x=sdate, y=proportion, shape=analysis, linetype=analysis)) + 
  geom_point() + 
  geom_smooth(method=lm, se=F, color="black") +
  facet_grid(.~ dv) + 
  scale_y_continuous(limits=c(0,1)) + 
  theme_light(16) + 
  labs(y = "Proportion Giving Response", 
       x = "Time",
       shape = "Response", 
       linetype = "Response") + 
  scale_x_date(date_breaks= "4 year", date_labels="%Y") + 
  theme(legend.position="top",
        legend.direction="horizontal", 
        legend.background = element_rect()) + 
  scale_shape_manual(values=c(17, 1, 0)) + 
  scale_linetype_manual(values=c(1, 3, 2)) + 
  theme(legend.position=c(0.30,0.95),
        legend.direction="horizontal", 
        legend.background = element_rect(), 
        legend.box.background = element_rect(colour = "black")) 

f2a <- f2 + 
  theme(legend.position = c(0.15, 0.9)) + 
   guides(shape = guide_legend(nrow=2))

 
ggsave(plot = f2a, "figure_oe2.png", 
       height=6, width=10, dpi=1200)


#########################################################################
                  # Figure OE3: Replication of Figure 2
#########################################################################

#####################################
#Left Hand Figure
#####################################

########Get weighted means
ideol_means <- comb_design %>%
  mutate(textreme_nodk = factor(textreme_nodk, 
                                levels=c(1,2,3), 
                                labels=c("Neither", "One Party", "Both Parties")), 
         ideol_ext = factor(ideol_ext, 
                            levels=c(1,2,3), 
                            labels=c("Moderate", "Liberal/Conservative", 
                                     "Very Liberal/Conservative"))) %>%
  filter(!is.na(textreme_nodk), !is.na(ideol_ext)) %>%
  group_by(sdate, ideol_ext, textreme_nodk) %>%
  summarize(proportion = srvyr::survey_mean(na.rm=T)) %>%
  mutate(lower = proportion - (1.96*proportion_se), 
         upper = proportion + (1.96*proportion_se), 
         sdate = as.Date(sdate))

#######Figure
f3 <- ideol_means %>%
  ggplot(aes(x=sdate, y=proportion, shape=textreme_nodk, linetype=textreme_nodk)) + 
  geom_point() + 
  geom_smooth(method=lm, se=F, color="black") + 
  facet_wrap(~ideol_ext) + 
  scale_y_continuous(limits=c(0,1)) + 
  theme_light(14) + 
  labs(title = "Mean Perceptions by Ideology (No DKs)", 
       y = "Proportion Giving Response", 
       x = "Time",
       shape = "Response", 
       linetype = "Response") + 
  scale_x_date(date_breaks= "4 year", date_labels="%Y") + 
  theme(legend.position = c(0.1, 0.9), 
        strip.text.x = element_text(face="bold", color="black"))  + 
  scale_shape_manual(values=c(17, 1, 0)) + 
  scale_linetype_manual(values=c(1, 3, 2))

f3a <- lemon::reposition_legend(f3, 'top left', offset= 0.005, 
                                panel = "panel-1-1")

#######################
#Mlogit Results
#From Stata
#######################

#Load, Clean, Bind the Data
neither_noc <- import("./Data/probs_1_nocontrol_nodk.dta") %>%
  rename(ame = `_margin`, 
         se = `_se`, 
         lower = `_ci_lb`, 
         upper = `_ci_ub`, 
         ideol = `_by1`) %>%
  mutate(outc = "Neither Party") %>%
  select(ame, se, lower, upper, ideol, outc)

one_noc <- import("./Data/probs_2_nocontrol_nodk.dta") %>%
  rename(ame = `_margin`, 
         se = `_se`, 
         lower = `_ci_lb`, 
         upper = `_ci_ub`, 
         ideol = `_by1`) %>%
  mutate(outc = "One Party") %>%
  select(ame, se, lower, upper, ideol, outc)

both_noc <- import("./Data/probs_3_nocontrol_nodk.dta") %>%
  rename(ame = `_margin`, 
         se = `_se`, 
         lower = `_ci_lb`, 
         upper = `_ci_ub`, 
         ideol = `_by1`) %>%
  mutate(outc = "Both Parties") %>%
  select(ame, se, lower, upper, ideol, outc)

ame_comb1 <- bind_rows(neither_noc, one_noc, both_noc) %>%
  mutate(outc = factor(outc, levels=c("Neither Party", "One Party", "Both Parties")), 
         ideol = factor(ideol, levels=c(1,2,3), 
                        labels=c("Moderate", "Liberal\n/Conservative", "Very Liberal\n/Conservative")))

#Figure
f3b <- ame_comb1 %>%
  ggplot(aes(y=ame, x=ideol, shape=outc)) + 
  geom_pointrange(aes(ymin=lower, ymax=upper)) + 
  facet_grid(. ~ outc ) + 
  theme_light(14) + 
  geom_hline(yintercept=0) + 
  labs(title = "Effect of Year by Ideology", 
       y = NULL, x = NULL) + 
  theme(legend.position = "none", 
        strip.text.x = element_text(face="bold", color="black"), 
        axis.text.x = element_text(size=7)) + 
  scale_shape_manual(values=c(17,1,0))

#Combine with f1a above and save
ggpubr::ggarrange(f3a, f3b , 
                  labels=c("A", "B"))

ggsave(".figure_oe3.png", 
       height=8, width=14, 
       dpi=1200)



#########################################################################
              # Figure OE4: Replication of Figure 3
#########################################################################


#############
#Pew/CNN
#############

#Import Data
pew_neither <- import("./Data/probs_1_nodk.dta")  %>%
  mutate(outcome = "Neither Party")

pew_one <- import("./Data/probs_2_nodk.dta")  %>%
  mutate(outcome = "One Party")

pew_both <- import("./Data/probs_3_nodk.dta")  %>%
  mutate(outcome = "Both Parties")

##Combine
pew_comb <- rbind(pew_neither, pew_one, pew_both) %>%
  rename(probs = `_margin`, 
         lower = `_ci_lb`, 
         upper = `_ci_ub`, 
         term = `_term`) %>%
  mutate(term = factor(term, 
                       levels=c(1,2), 
                       labels=c("PID Extremity", 
                                "Ideological Extremity"))) %>%
  select(-c(`_deriv`, `_predict`, `_at`, 
            `_atopt`, `_pvalue`, `_se`, `_statistic`))

##Figures

#Ideological Extremity
f41 <- pew_comb %>%
  filter(term == "Ideological Extremity") %>%
  rename(ideolext = `_m2`) %>%
  mutate(ideolext = factor(ideolext, levels=c(1,2,3), 
                           labels=c("Moderate", "Libs/Cons", "Very Libs/Cons")), 
         outcome = factor(outcome, levels=c("Neither Party", "One Party", "Both Parties"))) %>%
  ggplot(aes(x=ideolext, y=probs, group = outcome, shape = outcome)) + 
  geom_pointrange(aes(ymin=lower, ymax=upper)) + 
  geom_line() + 
  theme_bw() + 
  ylim(0,1) + 
  labs(title = "Pew/CNN: Ideological Extremity", 
       y = NULL, 
       x = NULL, 
       group = "Response", 
       shape = "Response", 
       x = "Ideological Extremity") + 
  theme(legend.position = "bottom", 
        legend.title=element_text(size=16), 
        legend.text=element_text(size=16)) + 
  ggpubr::theme_pubclean() + 
  scale_shape_manual(values=c(17, 1, 0)) 

#PID Extremity
f42 <- pew_comb %>%
  filter(term == "PID Extremity") %>%
  rename(pidext = `_m1`) %>%
  mutate(pidext = factor(pidext, levels=c(1,2,3), 
                         labels=c("Independent", "Leaning Partisan", "Partisan")), 
         outcome = factor(outcome, levels=c("Neither Party", "One Party", "Both Parties"))) %>%
  ggplot(aes(x=pidext, y=probs, group = outcome, shape = outcome)) + 
  geom_pointrange(aes(ymin=lower, ymax=upper)) + 
  geom_line() + 
  theme_bw() + 
  ylim(0,1) + 
  labs(title = "Pew/CNN: Party ID Extremity", 
       y = NULL, 
       x = NULL,
       group = "Response", 
       shape = "Response", 
       x = "PID Extremity") + 
  theme(legend.position = "bottom", 
        legend.title=element_text(size=16), 
        legend.text=element_text(size=16)) + 
  ggpubr::theme_pubclean() + 
  scale_shape_manual(values=c(17, 1, 0)) 


#############
#Lucid
#############
###Ideological Extremity
f43 <- import("./Data/lucid_sym_nodk.dta") %>%
  rename(outcome = `_predict`, 
         probs = `_margin`, 
         se = `_se`, 
         lower = `_ci_lb`, 
         upper = `_ci_ub`, 
         ideolext = `_m1`) %>%
  mutate(ideolext = factor(ideolext, 
                           levels=c(1,2,3,4), 
                           labels=c("Moderate.", "2", "3", "Very Lib/Cons")), 
         outcome = factor(outcome, 
                          levels=c(1,2,3), 
                          labels=c("Neither Party", 
                                   "One Party", 
                                   "Both Parties"))) %>%
  ggplot(aes(x=ideolext, y=probs, group=outcome, line=outcome, shape=outcome)) + 
  geom_pointrange(aes(ymin=lower, ymax=upper)) + 
  geom_line() + 
  theme_bw() + 
  labs(title = "Lucid: Ideological Extremity", 
       x = NULL,
       y = NULL, 
       shape = "Response", 
       line = "Response") + 
  scale_y_continuous(limits=c(0,1)) + 
  ggpubr::theme_pubclean() + 
  scale_shape_manual(values=c(17, 1, 0)) + 
  theme(legend.title=element_text(size=16), 
        legend.text=element_text(size=16))

####Operational Ideology
f44 <-import("./Data/lucid_op_nodk.dta") %>%
  rename(outcome = `_predict`, 
         probs = `_margin`, 
         se = `_se`, 
         lower = `_ci_lb`, 
         upper = `_ci_ub`, 
         opideolcons = `_at3`) %>%
  mutate(outcome = factor(outcome, 
                          levels=c(1,2,3), 
                          labels=c("Neither Party", 
                                   "One Party", 
                                   "Both Parties"))) %>%
  ggplot(aes(x=opideolcons, y=probs, group=outcome, line=outcome, shape=outcome)) + 
  geom_pointrange(aes(ymin=lower, ymax=upper)) + 
  geom_line() + 
  theme_bw() + 
  theme(legend.position = "bottom") + 
  labs(title = "Lucid: Issue Consistency", 
       x = NULL,
       y = NULL, 
       shape = "Response", 
       line = "Response") + 
  scale_y_continuous(limits=c(0,1)) + 
  ggpubr::theme_pubclean() + 
  scale_shape_manual(values=c(17, 1, 0)) + 
  theme(legend.title=element_text(size=16), 
        legend.text=element_text(size=16))


####PPD
f45 <-import("./Data/lucid_ppd_nodk.dta") %>% 
  rename(outcome = `_predict`, 
         probs = `_margin`, 
         se = `_se`, 
         lower = `_ci_lb`, 
         upper = `_ci_ub`, 
         ppd = `_at4`) %>%
  mutate(outcome = factor(outcome, 
                          levels=c(1,2,3), 
                          labels=c("Neither Party", 
                                   "One Party", 
                                   "Both Parties"))) %>%
  ggplot(aes(x=ppd, y=probs, group=outcome, line=outcome, shape=outcome)) + 
  geom_pointrange(aes(ymin=lower, ymax=upper)) + 
  geom_line() + 
  theme_bw() + 
  theme(legend.position = "bottom") + 
  labs(title = "Lucid: Perceived Policy Distance", 
       x = NULL,
       y = NULL, 
       shape = "Response", 
       line = "Response") + 
  scale_y_continuous(limits=c(0,1)) + 
  ggpubr::theme_pubclean() +
  scale_shape_manual(values=c(17, 1, 0)) + 
  theme(legend.title=element_text(size=16), 
        legend.text=element_text(size=16))

###PID
f46 <-import("./Data/lucid_pid_nodk.dta") %>%
  rename(outcome = `_predict`, 
         probs = `_margin`, 
         se = `_se`, 
         lower = `_ci_lb`, 
         upper = `_ci_ub`, 
         pididstr = `_at1`) %>%
  mutate(outcome = factor(outcome, 
                          levels=c(1,2,3), 
                          labels=c("Neither Party", 
                                   "One Party", 
                                   "Both Parties"))) %>%
  ggplot(aes(x=pididstr, y=probs, group=outcome, line=outcome, shape=outcome)) + 
  geom_pointrange(aes(ymin=lower, ymax=upper)) + 
  geom_line() + 
  theme_bw() + 
  theme(legend.position = "bottom") + 
  labs(title = "Lucid: Party Identity Strength", 
       x = NULL,
       y = NULL, 
       shape = "Response", 
       line = "Response") + 
  scale_y_continuous(limits=c(0,1)) + 
  scale_x_continuous(breaks=c(1,2,3,4), 
                     labels=c("Pure Ind.", "2", "3", "Strongest\nPartisan")) + 
  ggpubr::theme_pubclean() + 
  scale_shape_manual(values=c(17, 1, 0)) + 
  theme(legend.title=element_text(size=16), 
        legend.text=element_text(size=16))

#############
#Combined Figure
#############
ggpubr::ggarrange(f41, f43, f44, f45, f42, f46, 
                  common.legend=T, legend = "top", 
                  labels=c("A", "B", "C", "D", "E", "F"), 
                  font.label = list(size=18)) 


ggsave("figure_oe4.png", 
       height=8, width=14)

#########################################################################
          # Figure OD5: Replication of Figure 4
#########################################################################

###############
#Pew
###############

#Load Data
pew1  <- import("./Data/pewcnn_probs_alt_nodk.dta") %>%
  rename(probs = `_margin`, 
         se = `_se`, 
         zstat = `_statistic`, 
         pval = `_pvalue`, 
         lower = `_ci_lb`, 
         upper = `_ci_ub`, 
         ideol = `_at2`, 
         partisan = `_at1`, 
         outcome = `_predict`) %>%
  mutate(partisan = factor(partisan, 
                           levels=c(1,3), 
                           labels=c("Independent", "(Strong) Partisan")),
         ideol = factor(ideol, 
                        levels=c(1,2,3), 
                        labels=c("Moderate", "2", "Very Lib\nCons")), 
         outcome = factor(outcome, 
                          levels=c(1,2,3), 
                          labels=c("Neither Party", "One Party", "Both Parties"))) %>%
  select(probs, se, zstat, pval, lower, upper, ideol, partisan, outcome)


pew1a <- pew1 %>%
  ggplot(aes(x=ideol, y=probs, group=partisan, shape=partisan, linetype=partisan)) + 
  geom_pointrange(aes(ymin=lower, ymax=upper)) + 
  geom_line() + 
  facet_grid(.~outcome) + 
  labs(title = "Pew: Ideological Extremity", 
       y = NULL, x = NULL, shape = "Partisan", linetype="Partisan") +
  theme_pubclean() + 
  scale_y_continuous(limits=c(0,1)) + 
  scale_shape_manual(values=c(1,17))


###############
#Lucid
###############

#####Symbolic
luc_sym <- import("./Data/lucid_probs_ideolext_nodk.dta") %>%
  rename(probs = `_margin`, 
         se = `_se`, 
         zstat = `_statistic`, 
         pval = `_pvalue`, 
         lower = `_ci_lb`, 
         upper = `_ci_ub`, 
         ideol = `_at2`, 
         partisan = `_at1`, 
         outcome = `_predict`) %>%
  mutate(ideol = factor(ideol, 
                        levels=c(1,2,3,4), 
                        labels=c("Moderate" , "2", "3", "Very Lib\n/Cons")), 
         partisan=factor(partisan, 
                         levels=c(1,4), 
                         labels=c("Inependent", "(Strong) Partisan")), 
         outcome = factor(outcome, 
                          levels=c(1,2,3), 
                          labels=c("Neither Party", "One Party", "Both Parties"))) %>%
  select(probs, se, zstat, pval, lower, upper, ideol, partisan, outcome)

lucid_symb_figure <- luc_sym %>%
  ggplot(aes(x=ideol, y=probs, group=partisan, shape=partisan, linetype=partisan)) + 
  geom_pointrange(aes(ymin=lower, ymax=upper)) + 
  geom_line() + 
  facet_grid(.~outcome) + 
  labs(title = "Lucid: Ideological Extremity", 
       y = NULL, x = NULL, shape = "Partisan", linetype="Partisan") +
  theme_pubclean() + 
  scale_shape_manual(values=c(1,17)) + 
  scale_y_continuous(limits=c(0,1)) 

#####Operational
lucid_op <- import("./Data/lucid_probs_op_ideol_cons_nodk.dta") %>%
  rename(probs = `_margin`, 
         se = `_se`, 
         zstat = `_statistic`, 
         pval = `_pvalue`, 
         lower = `_ci_lb`, 
         upper = `_ci_ub`, 
         ideol = `_at2`, 
         partisan = `_at1`, 
         outcome = `_predict`) %>%
  mutate(ideol = factor(ideol, 
                        levels=c(0,1,2,3,4,5), 
                        labels=c("Mix.", "1", "2", "3", "4", "Cons.")), 
         partisan=factor(partisan, 
                         levels=c(1,4), 
                         labels=c("Inependent", "(Strong) Partisan")), 
         outcome = factor(outcome, 
                          levels=c(1,2,3), 
                          labels=c("Neither Party", "One Party", "Both Parties"))) %>%
  select(probs, se, zstat, pval, lower, upper, ideol, partisan, outcome)

lucid_op_fig <- lucid_op %>%
  ggplot(aes(x=ideol, y=probs, group=partisan, shape=partisan, linetype=partisan)) + 
  geom_pointrange(aes(ymin=lower, ymax=upper)) + 
  geom_line() + 
  facet_grid(.~outcome) + 
  labs(title = "Lucid: Operational Consistency", 
       y = NULL, x = NULL, shape = "Partisan", linetype="Partisan") +
  theme_pubclean() + 
  scale_shape_manual(values=c(1,17)) + 
  scale_y_continuous(limits=c(0,1))


####PPD
lucid_ppd <- import("./Data/lucid_probs_ppd_nodk.dta") %>%
  rename(probs = `_margin`, 
         se = `_se`, 
         zstat = `_statistic`, 
         pval = `_pvalue`, 
         lower = `_ci_lb`, 
         upper = `_ci_ub`, 
         ideol = `_at2`, 
         partisan = `_at1`, 
         outcome = `_predict`) %>%
  mutate(ideol = factor(ideol, 
                        levels=c(0,1,2,3,4,5,6), 
                        labels=c("Equal\nProx", "1", "2", "3", "4", "5", "Max R.\nProx" )), 
         partisan=factor(partisan, 
                         levels=c(1,4), 
                         labels=c("Inependent", "(Strong) Partisan")), 
         outcome = factor(outcome, 
                          levels=c(1,2,3), 
                          labels=c("Neither Party", "One Party", "Both Parties"))) %>%
  select(probs, se, zstat, pval, lower, upper, ideol, partisan, outcome)

lucid_ppd_fig <- lucid_ppd %>%
  ggplot(aes(x=ideol, y=probs, group=partisan, shape=partisan, linetype=partisan)) + 
  geom_pointrange(aes(ymin=lower, ymax=upper)) + 
  geom_line() + 
  facet_grid(.~outcome) + 
  labs(title = "Lucid: Perceived Policy Distance", 
       y = NULL, x = NULL, shape = "Partisan", linetype="Partisan") +
  theme_pubclean() + 
  scale_shape_manual(values=c(1,17)) + 
  scale_y_continuous(limits=c(0,1)) + 
  scale_x_discrete(breaks=c("Equal\nProx", "Max R.\nProx"))


###############
#Combining
###############
ggarrange(pew1a, lucid_symb_figure, 
          lucid_op_fig, 
          lucid_ppd_fig, 
          labels=c("A", "B", "C", "D"), 
          common.legend = T)


ggsave("figure_oe5.png", 
       height=8, width=12, dpi=1200)

